home *** CD-ROM | disk | FTP | other *** search
/ Almathera Ten Pack 2: CDPD 1 / Almathera Ten on Ten - Disc 2: CDPD 1.iso / pd / 076-100 / 085 / plot6 / plot6.c < prev    next >
C/C++ Source or Header  |  1995-03-13  |  6KB  |  239 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3.  
  4. struct ent{
  5.     float ra,dc;
  6.     short mg,cl;
  7. };
  8.  
  9. extern int MAX_X,MAX_Y;
  10.  
  11. /* 18.5cm/screen Y 26.2cm/screen X -> aspect ratio 1.4162 */
  12.  
  13. #define ASPECT 1.4162
  14.  
  15. /* MATH!
  16.  
  17.    spherical to rectangular coordinates.
  18.    earth's axis is z axis
  19.    dc measured from positive (north) z axis
  20.    ra measured from x axis (we'll fix that later)
  21.    All stars are on a unit sphere.
  22.    
  23.    x = sin(pi/2-dc)cos(ra)
  24.    y = sin(pi/2-dc)sin(ra)
  25.    z = cos(pi/2-dc)
  26.    
  27.    rotation rule (about origin [ie about z axis])
  28.    
  29.    a' = a cos(th) - b sin(th)
  30.    b' = a sin(th) + b cos(th)
  31.    
  32.    rotate to see desired line of longitude(ra0).
  33.    ie line up y axis with chosen longitude (remember we were measuring from x)
  34.    rotate counter clockwise pi/2, rotate clockwise ra0 about z (ie pi/2-ra0)
  35.  
  36.    x' = sin(pi/2-dc)cos(ra)cos(pi/2-ra0) - sin(pi/2-dc)sin(ra)sin(pi/2-ra0)
  37.       = sin(pi/2-dc)cos(ra+pi/2-ra0) [believe me!]
  38.    y' = sin(pi/2-dc)cos(ra)sin(pi/2-ra0) + sin(pi/2-dc)sin(ra)cos(pi/2-ra0)
  39.       = sin(pi/2-dc)sin(ra+pi/2-ra0)
  40.    z' = cos(pi/2-dc)
  41.    
  42.    now rotate down from the north pole along your longitude to desired latitude
  43.    ie rotate down pi/2 then up by dc0
  44.    we are rotating around the x axis
  45.    
  46.    x" = sin(pi/2-dc)cos(ra+pi/2-ra0)
  47.    y" = sin(pi/2-dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-cos(pi/2-dc)sin(pi/2-dc0)
  48.       = cos(dc)sin(ra+pi/2-ra0)cos(pi/2-dc0)-sin(dc)sin(pi/2-dc0)
  49.    z" = sin(pi/2-dc)sin(ra+pi/2-ra0)sin(pi/2-dc0)+cos(pi/2-dc)cos(pi/2-dc0)
  50.       = cos(dc)sin(ra+pi/2-ra0)sin(pi/2-dc0) + sin(dc)cos(pi/2-dc0)
  51.    
  52.    I set dc0 to pi/2-dc0 and ra0 ot pi/2-ra0.
  53.    a = sin(dc)        b = cos(dc)
  54.    c = sin(dc0)        d = cos(dc0)
  55.    e = sin(ra + ra0)    f = cos(ra + ra0)
  56.    
  57.    x" = df;
  58.    y" = deb-ca;
  59.    z" = dea+cb;
  60.    
  61.    Now! viewing transforms.  I assume we are at (0,0,0) looking at ra0 x dc0.
  62.    We are looking down the z axis with y up.
  63.    
  64.    our famous Ys = d/s Ye/Ze view box transform gives us:
  65.    (ie Y in screen coordinates from eye coordinates)
  66.    
  67.    Xs = df/(dea+cb);
  68.    Ys = (deb-ca)/(dea+cb);
  69.    
  70.    first we discard any stars behind us (ie z < 0)
  71.    if the star is more than half a screen to the side, discard it.
  72.    NB: we need not calculate f until the last moment.
  73.    
  74.    divide x by the aspect ratio, then map to physical coordinates.
  75.    
  76.    the rest is a subjective star plotting mechanism
  77. */
  78.  
  79. #define mabs(x) (((x)<0.0)?-(x):(x))
  80. #define mfloor(x) (((x)<0.0)?-(int)-(x):(int)(x))
  81.  
  82.  
  83. main(argc,argv)
  84. int argc;
  85. char **argv;
  86. {
  87.     FILE *fp;
  88.     struct ent *ar,*star;
  89.     int t,i,num;
  90.     double a,b,c,d,e,f,de,ra0,de0,wi,xt,yt;
  91.     double HLF;
  92.     double z;
  93.     double TO_RADS;
  94.     
  95.     if (argc <= 1) {
  96.         printf("usage plot <file-list>\n");
  97.     exit(0);
  98.     }
  99.     
  100.     /* beware the ffp constant expression bug! */
  101.     TO_RADS  = (TO_RADS=2.0) * 3.1415926535 / (TO_RADS=360.0);
  102.     
  103.     HLF = 90.0 * TO_RADS;
  104.     
  105.     printf("ra,de,wi [degrees]:");fflush(stdout);
  106.     scanf("%lf %lf %lf",&ra0,&de0,&wi);
  107.     ra0 = HLF-ra0*TO_RADS;
  108.     de0 = HLF-de0*TO_RADS;
  109.     wi = tan(mabs(wi)*TO_RADS);
  110.  
  111.     window_init();
  112.     a = sin(de0); b = cos(de0);
  113.  
  114.     while (argc-- > 1){
  115.         fp = fopen(argv[argc],"r");
  116.     
  117.     if (fp == NULL){
  118.         printf("cant find %s\n",argv[argc]);
  119.         exit(0);
  120.     }
  121.     
  122.     t = fread((char*)&num,4,1,fp);
  123.     if (t != 1) {
  124.         printf("didn't read file size!\n");
  125.         exit(0);
  126.     }
  127.     
  128.     printf("reading %d records\n",num);
  129.     fflush(stdout);
  130.     
  131.     ar = (struct ent*)malloc(num*sizeof(struct ent));
  132.     if (ar == NULL){
  133.         printf("couldn't malloc %d bytes\n",num*sizeof(struct ent));
  134.         exit(0);
  135.     }
  136.     
  137.     t = fread((char*)ar,sizeof(struct ent),num,fp);
  138.     fclose(fp);
  139.     
  140.     printf("read %d\n",t);
  141.     
  142.     for (i = 0;i<t;i++) {
  143.         check_input();
  144.         star = ar+i;
  145.         
  146.         c = sin(star->dc);
  147.         d = cos(star->dc);
  148.         
  149.         e = sin(star->ra + ra0);
  150.         
  151.         de = d * e;
  152.         z = wi*(de * a + c * b);
  153.         if (z <= 1.0e-16) continue;    /* behind us, so skip it */
  154.         
  155.         yt = (de * b - c * a) / z;
  156.         if (yt < -0.5 || yt > 0.5) continue;
  157.         
  158.         f = cos(star->ra + ra0);
  159.         
  160.         xt = d * f / ASPECT / z;
  161.         if (xt < -0.5 || xt > 0.5) continue;
  162.         
  163.         plot_star((int)((xt + 0.5)*MAX_X),
  164.               (int)((yt + 0.5)*MAX_Y),star->mg,star->cl);
  165.     }
  166.     printf("done\n");
  167.     fflush(stdout);
  168.     free (ar);
  169.     
  170.     }
  171.     while (1){
  172.         wt();
  173.     check_input();
  174.     }
  175. }
  176.  
  177. plot_star(x,y,mg,cl)
  178. int x,y;
  179. short mg,cl;
  180. {
  181.     int c,colour,m;
  182.     
  183.     m = (int)mfloor(((double) mg)/100.0+0.49);
  184.     
  185.     c = (int)mfloor(cl/60.0+0.49) + 1;
  186.  
  187.     if (c < 0) c = 0;
  188.     else if (c > 3) c = 3;
  189.     
  190.     c = 1 + 3 * c + 2;/*dimmest of range*/
  191.     
  192.     colour = c-1;/* medium */
  193.     
  194.     switch(m){
  195.         case  0: point(x-3,y,colour);
  196.              point(x-4,y,colour);
  197.          point(x+3,y,colour);
  198.          point(x+4,y,colour);
  199.          point(x+5,y,colour);
  200.          point(x-5,y,colour);
  201.          colour = c-2;
  202.          point(x+2,y+1,colour);
  203.          point(x+2,y-1,colour);
  204.          point(x-2,y+1,colour);
  205.          point(x-2,y-1,colour);
  206.     case 1:     point(x,y+2,colour);
  207.          point(x,y-2,colour);
  208.          colour = c-2;
  209.     case  2: /* in bright colour */
  210.          point(x+3,y,colour);
  211.          point(x-3,y,colour);
  212.          point(x+2,y,colour);
  213.          point(x-2,y,colour);
  214.          colour = c-2;
  215.          point(x+1,y+1,colour);
  216.          point(x+1,y-1,colour);
  217.          point(x-1,y+1,colour);
  218.          point(x-1,y-1,colour);
  219.          /* in brightest colour */
  220.     case  3: /* in dimmest colour */
  221.          colour = c-2;
  222.          /* brighten colour */
  223.     case  4: /* in dimmest colour*/
  224.          point(x,y+1,colour);
  225.          point(x,y-1,colour);
  226.          point(x+1,y,colour);
  227.              point(x-1,y,colour);
  228.          point(x+2,y,colour);
  229.              point(x-2,y,colour);
  230.     case  5: colour = c-2;/* set to brightest */
  231.          break;
  232.     case  6: colour = c-1;
  233.          break;
  234.         default: colour = c;/* in dimmest colour */
  235.     }
  236.     point(x,y,colour);
  237.  
  238. }
  239.